Down stream analysis of the dataset ‘1_IMR90_domains_23plex’
library(ggforce)
library(tidyverse)
dcolor <- ggsci::scale_color_d3('category20')
dfill <- ggsci::scale_fill_d3('category20')
ccolor <- viridis::scale_color_viridis()
style <- function(x,color,method='UMAP',axis=c(1,2),Theme=NULL,
axis.text=FALSE,axis.title=FALSE,coord_fix =TRUE,
palette=TRUE,title=NULL,
legend=TRUE,scale_color_log10=FALSE,
guide_pointSize=2.5,...) {
axis <- paste0(method,axis)
if(is.null(Theme)) Theme <- theme_bw
g <- ggplot() + Theme()
g <- g + geom_point(aes_(x=as.name(axis[1]),y=as.name(axis[2]),color=as.name(color)),x,...)
if(coord_fix) {
g <- g + coord_fixed()
}
if(!axis.text) {
g <- g + theme(axis.text = element_blank(),axis.ticks = element_blank())
}
if(!axis.title) {
g <- g + theme(axis.title = element_blank())
}
if(is.null(title)) {
g <- g + ggtitle(method)
}else if(!is.na(title)) {
g <- g + ggtitle(title)
}
if(!legend) g <- g + theme(legend.position = 'none')
if(scale_type(as.matrix(x[,color]))=='discrete' & !is.null(guide_pointSize)) {
g <- g + guides(color = guide_legend(override.aes = list(size=guide_pointSize)))
}
if(palette) {
if(scale_type(as.matrix(x[,color]))=='discrete') {
g <- g + ggsci::scale_color_d3('category20')
}else{
if(scale_color_log10){
g <- g + viridis::scale_color_viridis(trans='log10')
}else{
g <- g + viridis::scale_color_viridis()
}
}
}
return(g)
}
tib2df <- function(tib,RowNames=1) {
rn <- pull(tib,all_of(RowNames))
df <- tib %>%
select(-RowNames) %>%
as.data.frame()
rownames(df) <- rn
return(df)
}
data <- tibble(
path = list.files('1_IMR90_domains_23plex/Raw/',full.names = T),
sampleName = paste(rep(c('G','S'),each=5),formatC(1:10,width = 2,flag = 0),sep = '-')
) %>%
separate(sampleName,into=c('type','cellID'),remove = F) %>%
mutate(metric = map(path,read_csv,show_col_types=F)) %>%
select(-path) %>% unnest(metric) %>%
unite(Domain,sampleName,domainID,sep = '-',remove = F)
tmp <- data %>%
gather(key=feature,value=exp,-(1:5)) %>%
separate(feature,into=c("region","ab"),sep = " ") %>%
mutate(region = factor(region,levels=c("In","Mid","Ex")),) %>%
split(f=.$ab)
map(names(tmp),~{
tmp[[.x]] %>%
ggplot(aes(region,exp,color=cellID)) +
theme_bw() + geom_violin(color="grey") + geom_jitter(size=.6) +
facet_wrap(~type) + ggtitle(.x) + dcolor +
labs(y = 'Expression Level')
}) %>%
patchwork::wrap_plots(guides = "collect",ncol = 6)
X <- as.matrix(tib2df(data[,-(2:5)]))
p <- prcomp(X)
screeplot(p,type='l',n=30)
total <- rowSums(data[,-(1:5)])
cor(p$x[,1:30],total) %>%
plot(xlab = 'PC',ylab = 'CorCoef_vsTotalExp')
set.seed(10)
um <- uwot::umap(p$x[,2:10],metric = "cosine",pca = 10,n_neighbors = 10,min_dist = .3)
Meta <- data[,1:5] %>%
mutate(UMAP1=um[,1],UMAP2=um[,2])
style(Meta,"type",size=1,palette=F)
style(Meta,"cellID",size=1)
set.seed(1)
D <- as.matrix(dist(p$x[,2:10]))
k <- 10
adj <- apply(D,1,function(x) rank(x)%in%2:(k+1))
g <- igraph::graph_from_adjacency_matrix(adj,mode='undirected')
lei <- igraph::cluster_leiden(g,resolution_parameter = .02)
Meta <- Meta %>% mutate(cluster=as.factor(lei$membership))
style(Meta,"cluster") + ggsci::scale_color_lancet()
Scale for colour is already present.
Adding another scale for colour, which will replace the existing scale.
tmp <- data %>% select(-sampleName,-type,-cellID,-domainID) %>%
tib2df() %>% scale %>% t %>% scale %>% asinh()
idx <- rownames(tmp) %>% sub(".* ","",.)
mod <- model.matrix(~0+idx)
colnames(mod) <- sub('idx','',colnames(mod))
tmp2 <- (t(tmp)%*%mod)/3
Dists <- dist(t(tmp2))
h <- hclust(Dists, method = "average")
dend <- as.dendrogram(h)
ggdend <- ggdendro::dendro_data(dend)
ggdend$segments$y <- ggdend$segments$y/3
ggdend$segments$yend <- ggdend$segments$yend/3
g <- ggplot(ggdendro::segment(ggdend))
g + geom_segment(aes(x=y,y=x,xend=yend,yend=xend)) + scale_x_reverse() +
geom_text(data = ggdend$label, aes(y, x, label = label), hjust = 0, angle = 0) + xlim(10,-8) + theme_void()
Scale for x is already present.
Adding another scale for x, which will replace the existing scale.
n <- ggdend$labels$label %>% length
idx <- ggdend$labels$label %>% rev %>% rep(each=3) %>% paste0(rep(c("In ","Mid ","Ex "),n),.)
tmp <- tmp[idx,]
anno <- data %>% select(Domain,sampleName) %>%
separate(sampleName,into = c('type','cell')) %>%
tib2df()
anno$cluster <- as.factor(lei$membership)
anno.color <- list()
anno.color$type <- c('G' = "#F8766D", 'S' = "#00BFC4")
anno.color$cell <- ggsci::scale_color_d3()$palette(10)
anno.color$cluster <- ggsci::scale_color_lancet()$palette(4)
names(anno.color$cell) <- formatC(1:10,width = 2,flag = '0')
names(anno.color$cluster) <- 1:4
pheatmap::pheatmap(tmp[,lei$membership == 1],
cluster_rows = F,
annotation_col = anno,annotation_colors = anno.color,show_colnames = F,
color = viridis::viridis(32),border_color = NA,cellheight = 10,cellwidth = 2.5)
pheatmap::pheatmap(tmp[,lei$membership == 2],
cluster_rows = F,
annotation_col = anno,annotation_colors = anno.color,show_colnames = F,
color = viridis::viridis(32),border_color = NA,cellheight = 10,cellwidth = 2.5)
pheatmap::pheatmap(tmp[,lei$membership == 3],
cluster_rows = F,
annotation_col = anno,annotation_colors = anno.color,show_colnames = F,
color = viridis::viridis(32),border_color = NA,cellheight = 10,cellwidth = 2.5)
pheatmap::pheatmap(tmp[,lei$membership == 4],
cluster_rows = F,
annotation_col = anno,annotation_colors = anno.color,show_colnames = F,
color = viridis::viridis(32),border_color = NA,cellheight = 10,cellwidth = 2.5)
sessionInfo()
R version 4.2.2 (2022-10-31)
Platform: aarch64-apple-darwin20 (64-bit)
Running under: macOS Ventura 13.2.1
Matrix products: default
LAPACK: /Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/lib/libRlapack.dylib
locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] SeuratObject_4.1.3 Seurat_4.3.0.1 ggforce_0.4.1 lubridate_1.9.2 forcats_1.0.0 stringr_1.5.0 dplyr_1.1.2 purrr_1.0.1 readr_2.1.4
[10] tidyr_1.3.0 tibble_3.2.1 ggplot2_3.4.2 tidyverse_2.0.0
loaded via a namespace (and not attached):
[1] utf8_1.2.3 spatstat.explore_3.2-1 reticulate_1.30 tidyselect_1.2.0 RSQLite_2.3.1 AnnotationDbi_1.60.2 htmlwidgets_1.6.2
[8] grid_4.2.2 BiocParallel_1.32.6 Rtsne_0.16 munsell_0.5.0 codetools_0.2-19 ragg_1.2.5 ica_1.0-3
[15] future_1.33.0 miniUI_0.1.1.1 withr_2.5.0 spatstat.random_3.1-5 colorspace_2.1-0 progressr_0.13.0 Biobase_2.58.0
[22] knitr_1.43 rstudioapi_0.15.0 stats4_4.2.2 ROCR_1.0-11 tensor_1.5 listenv_0.9.0 labeling_0.4.2
[29] GenomeInfoDbData_1.2.9 polyclip_1.10-4 bit64_4.0.5 farver_2.1.1 pheatmap_1.0.12 parallelly_1.36.0 vctrs_0.6.3
[36] generics_0.1.3 xfun_0.39 timechange_0.2.0 R6_2.5.1 phateR_1.0.7 GenomeInfoDb_1.34.9 bitops_1.0-7
[43] spatstat.utils_3.0-3 cachem_1.0.8 promises_1.2.0.1 scales_1.2.1 vroom_1.6.3 gtable_0.3.3 globals_0.16.2
[50] goftest_1.2-3 rlang_1.1.1 systemfonts_1.0.4 RcppRoll_0.3.0 pamr_1.56.1 splines_4.2.2 lazyeval_0.2.2
[57] spatstat.geom_3.2-4 yaml_2.3.7 reshape2_1.4.4 abind_1.4-5 stylo_0.7.4 httpuv_1.6.11 tools_4.2.2
[64] tcltk_4.2.2 ellipsis_0.3.2 jquerylib_0.1.4 RColorBrewer_1.1-3 ggdendro_0.1.23 proxy_0.4-27 BiocGenerics_0.44.0
[71] ggridges_0.5.4 Rcpp_1.0.11 plyr_1.8.8 zlibbioc_1.44.0 RCurl_1.98-1.12 deldir_1.0-9 pbapply_1.7-2
[78] viridis_0.6.4 cowplot_1.1.1 S4Vectors_0.36.2 zoo_1.8-12 ggrepel_0.9.3 cluster_2.1.4 magrittr_2.0.3
[85] data.table_1.14.8 scattermore_1.2 lmtest_0.9-40 RANN_2.6.1 fitdistrplus_1.1-11 Signac_1.10.0 matrixStats_1.0.0
[92] hms_1.1.3 patchwork_1.1.2 mime_0.12 evaluate_0.21 xtable_1.8-4 IRanges_2.32.0 gridExtra_2.3
[99] compiler_4.2.2 KernSmooth_2.23-22 crayon_1.5.2 htmltools_0.5.5 later_1.3.1 tzdb_0.4.0 DBI_1.1.3
[106] tweenr_2.0.2 formatR_1.14 MASS_7.3-60 Matrix_1.5-4.1 cli_3.6.1 parallel_4.2.2 igraph_1.5.0.1
[113] GenomicRanges_1.50.2 pkgconfig_2.0.3 sp_2.0-0 plotly_4.10.2 spatstat.sparse_3.0-2 bslib_0.5.0 XVector_0.38.0
[120] digest_0.6.33 sctransform_0.3.5 RcppAnnoy_0.0.21 tsne_0.1-3.1 spatstat.data_3.0-1 Biostrings_2.66.0 rmarkdown_2.23
[127] leiden_0.4.3 fastmatch_1.1-3 uwot_0.1.16 shiny_1.7.4.1 Rsamtools_2.14.0 lifecycle_1.0.3 nlme_3.1-162
[134] jsonlite_1.8.7 viridisLite_0.4.2 fansi_1.0.4 pillar_1.9.0 ggsci_3.0.0 lattice_0.21-8 KEGGREST_1.38.0
[141] fastmap_1.1.1 httr_1.4.6 survival_3.5-5 glue_1.6.2 png_0.1-8 bit_4.0.5 tcltk2_1.2-11
[148] class_7.3-22 stringi_1.7.12 sass_0.4.7 blob_1.2.4 textshaping_0.3.6 memoise_2.0.1 irlba_2.3.5.1
[155] e1071_1.7-13 future.apply_1.11.0 ape_5.7-1